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ABSTRACT 

We model for the first time the complete orbital evolution of a pair of Supermassive Black Holes 
(SMBHs) in a 1:10 galaxy merger of two disk dominated gas-rich galaxies, from the stage prior to the 
formation of the binary up to the onset of gravitational wave emission when the binary separation has 
shrunk to 1 milli parsec. The high-resolution smoothed particle hydrodynamics (SPH) simulations 
used for the first phase of the evolution include star formation, accretion onto the SMBHs as well 
as feedback from supernovae explosions and radiative heating from the SMBHs themselves. Using 
the direct A-body code 0-GPU we evolve the system further without including the effect of gas, 
which has been mostly consumed by star formation in the meantime. We start at the time when 
the separation between two SMBHs is ~ 700 pc and the two black holes are still embedded in their 
galaxy cusps. We use 3 million particles to study the formation and evolution of the SMBH binary 
till it becomes hard. After a hard binary is formed, we reduce (reselect) the particles to 1.15 million 
and follow the subsequent shrinking of the SMBH binary due to 3-body encounters with the stars. 
We find approximately constant hardening rates and that the SMBH binary rapidly develops a high 
eccentricity. Similar hardening rates and eccentricity values are reported in earlier studies of SMBH 
binary evolution in the merging of dissipationless spherical galaxy models. The estimated coalescence 
time is ~ 2.9 Gyr, significantly smaller than a Hubble time. We discuss why this timescale should be 
regarded as an upper limit. Since 1:10 mergers are among the most common interaction events for 
galaxies at all cosmic epochs, we argue that several SMBH binaries should be detected with currently 
planned space-borne gravitational wave interferometers, whose sensitivity will be especially high for 
SMBHs in the mass range considered here. 



Subject headings: black hole physics - galaxies: evolution 
gravitational waves. 



galaxies: interactions - galaxies: nuclei 



1. INTRODUCTION 

Central supermassive black holes (SMBHs) are ubiq- 
uitous and are found in a variety of galaxies, ranging 
from low mass galaxies to the most massive early-typ e 
galaxies (jFerrarese fc Fordl 12001 iGultekin et al.1 12009). 
Within our current cosmological picture of hierarchical 
structure formation, galaxies form through continuous 
mergers. If both candidate galaxies harbor a central 
SMBH before the merger, the evolution of the latter 
is thought to be as follows (Bcgclm an et al.lll980f ): the 
SMBHs of the merging galaxies sink towards the cen- 
ter of the merger remnant due to dynamical friction 
and form a gravitationally bound binary system. The 
further evolution of the SMBH binary is governed by 
interactions with stars and gas. If the binary semi- 
major axis value shrinks to a value where emission of 
gravitational waves (GWs) efficiently takes away energy 
and angular momentum from the binary, the coalescence 
of SMBHs becomes inevitable. In fact there is grow- 
ing observational evidence for this process: there are 
reports about t wo widely separated SMBHs in a sin- 
gle galaxy (e.g. iKomossa et. al. 1 120031 : [Rodrigue z et al.1 



2006; iFabbiano et al.l 120111 ) as well as indirect ev i- 
dence for bina r y SMBHs (e.g . [Merritt fc Ekers 200 
Liu et all [20031: IValtonen et al.l [20081: iBoroson fc Laue 



20091: llguchi'et al. 2010). 



Binary SMBHs are of particular interest in astro- 
physics and general relativity: if SMBH binary evo- 
lution leads to coalescence following a galaxy merger, 
such an event would give rise to one of the loud- 
est possible bursts of GWs detectable fo r future space 
borne low-frequency laser interferometers ([Hughes! [20031 : 
iBarack fc Cutlerll2004t[Amaro-Seoane et al.ll2012f )7 How- 
ever, the exact pathway to SMBH coalescence and espe- 
cially the associated timescales, and hence the chances 
for a possible detection during any satellite mission run- 
time, is still a matter of debate. 

Numerical A-body simulations considering the de- 
cay of a pair of SMBHs in major mergers of mas- 
sive elliptical galaxies due to the interaction with 
the stellar background only show that the critical 
separation for gravitational wave emission should be 
reached in about one Gyr as long as the galaxies de- 
viate sufficiently from sphericity, a configuration which 
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allow s to have centrophilic o r bits dSridhar fe Toumal 
l99l I Poon fc Merritti J 200lt IMerritt fc Poonl 120041: 
Merritt fc Vasilievl I2011J) that efficiently refi l l the 
loss cone (Bcrczi k et al.l 120 06: Bc rentzen et al.l 120091: 



Khanet all I2011L 12012: Pr eto et al.1 120111 ). iKhan et al.1 
(2012) showed that SMBHs with masses ranging from 
10 6 M© — 10 7 Mq, e.g., at the center of faint elliptical 
galaxies or in the bulge of spiral galaxies, coalesce well 
within a Gyr after the merger of the two galaxies/bulges. 
However, these simulations started from idealized initial 
conditions and do not consider the effects of gas dynam- 
ics, star formation, etc. 

Conversely, numerical simulations investigating major 
mergers including hydrodynamics show that a SMBH 
binary can form very rapidly after the galaxy colli- 
sion takes place, sometimes in less than a million years 
(jMaver et a l. 2007), but are not yet conclusive on the 
subsequent shrinking of the binary at sub-parsec scale 
separations due to the prohibitive cost of the compu- 
tation with increasing resolution and the uncertainties 
in the modeling of g as thermodynamics and turbulence 
(|Chapon et al.ll201l[) . 

Furthermore, minor mergers are the most frequent 
type of mergers in the Universe, especially those with 
a mass ratio of some 1 : 10. In addition, they usu- 
ally involve disk-dominated galaxies because these are 
not on ly the most commo n galaxies in the Universe 
today (Binggcli et al. 1988) but they were even more 
common at higher redshift since they are the likely 
progenitors of present-da y massive early type galaxies 
(IVan der Wel"^t~aT1[20rl . 

Simulations of minor mergers between disk galaxies 
show that the role of gas can be crucial in delivering 
the pair of SMBHs to the center of the merger rem- 
nant; without gas the SMBH of the secondary galaxy 
is left wandering at kilo-parsec distances, at which the 
dynamical friction timescale is longer than the Hubble 
time, because the core of the galaxy that is hosting it is 
tidally disrup ted before it can sink to the center of the 
primary (Callcgari et al.l [20091 ). With the gas the pair- 
ing of the two SMBHs is successful because the gas-rich 
merging satellite undergoes a strong central star forma- 
tion burst, developing a higher central density allowing 
its survi val to tidal disruption down to the center of the 
primary. iCallegari et al.1 (J2011I ) , who studied the pairing 
of SMBHs in the merger of late type gas rich galaxies 
in simulations which include star formation and accre- 
tion onto the SMBH, found that the mass ratio of the 
two SMBHs can change significantly due to the fact that 
the secondary is fed at a higher rate as a result of the 
stronger gas inflow caused by the tidal disturbance. The 
latter simulations of minor mergers could not follow the 
decay of the two SMBHs to separations of less than tens 
of parsec. Hence it is unknown whether or not they would 
reach the critical separation for gravitational wave emis- 
sion and on which timescales. Continuing the calcula- 
tions further with gas dynamics is computationally chal- 
lenging. In addition, at the end of the merger significant 
star formation takes place, increasing the contribution of 
the stars to the potential relative to the gas. 

Therefore it is sensible to consider the approximation 
in which the stars are the dominant source of the drag 
acting onto the SMBHs during the last stage of the sink- 




FlG. 1. — Density distribution of the dark matter (top panels) 
and stellar component (bottom panels) in the x — y (left column) 
and y — z (right column) planes. The two high density regions are 
clearly visible around the two SMBHs (black dots) in the center. 
The size of each box is 4 kpc. 

ing before gravitational waves take over. This is what 
we attempt to do for the first time in this paper, start- 
ing from a realistic la te configuration of a g as-rich minor 
merger performed bv ICallegari et al.1 (|2011l ). we followed 
the evolution of the SMBH binary to the milli parsec 
(mpc) regime where GWs start becoming important. 

The paper is organized in the following way: In Sec- 
tion 2 we describe the initial conditions for our direct 
iV-body simulations. The numerical methods and codes 
used for our simulations are described in Section 3. In 
Section 4 we explain the results for the evolution of the 
SMBH binary. Finally, Section 5 concludes the findings 
of our study. 

2. INITIAL CONDITIONS 

Owing to numerical limitations in simulations of galaxy 
mergers in respect of accuracy on parsec or sub-parsec 
scale, e.g., due to approximations in the force calculation 
(gravitational softening, tree-approximation) or insuffi- 
cient integration schemes we follow a multi-step strategy 
to study the evolution of binary SMBHs. We take sim- 
ulation data from a galaxy merger simulation including 
SMBHs in its final phase, i.e., before reaching resolu- 
tion limits and then follow up these simulations in high- 
resolution using direct ./V-body calculations with the dy- 
namically relevant though spatially truncated region of 
t he original galaxy merg er . 

ICallegari et al.l (|2011l ) studied the pairing of binary 
SMBHs in several 1:10 mergers of two disk galax- 
ies using iV-body/Smoothed Particle Hydrodynamics 
SPH) simulations co nducted with the GASOLINE code 
Wadsley et al.l 120041 ) . The reference galaxy model in 
that study was a Milky Way type disk galaxy consist- 
ing of three components: (i) a spherical and isotropic 
Navarro, Frenk & White profile dark matter halo, (ii) an 
exponential disk consisting of stars and gas and (iii) a 
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spherical Hernquist bulge. The mergers considered dif- 
ferent orbital configurations. In all cases, except for a pe- 
culiar retrograde merger case, the two SMBHs paired at 
the center of the remnant in less than 1 Gyr from the end 
of the galaxy collision. In this paper we employ one par- 
ticular merger set-up as initial conditions since already a 
single simulation of this kind requires several months of 
computation. In particular we choose the coplanar pro- 
grade merger simulation with 30% gas fraction in the 
disk since it resulted in the shortest merging time for the 
galaxies. 

The gravitat i onal s oftening adopted in the study of 
ICallegari et al.l ()201lD was e = 45 pc for dark matter 
and baryonic particles in the larger galaxy, whereas for 
smaller galaxy it was 20 pc. The SMBH was intro- 
duced at the center of each galaxy represented by a point 
mass particle. The adopted masses for the SMBHs were 
6 x 10 5 M Q and 6 x 10 4 M Q for the primary and satel- 
lite gala xy respectively, cons istent with the M m - Mb u io e 
relation (jMerritt fc Ferraresdl20"o"lt lHaring fc Rixll2004f) . 
The simulation also includes the effects of star formation 
and accretion on to the SMBHs, as well as feedback from 
both processes. The final separation of the two SMBHs 
at the end of the simulation was ~ 30 pc, which is com- 
parable to the softening that was used in the simulations. 
For more details of the construction of the galaxy models 
and the set up for the initial orbit of the g alaxy merger we 
refer the reader to Call egari et al.l (j2011h . Here it suffices 
to say that the merger was assumed to start at z ~ 3 — 5, 
and is completed at a time corresponding to z ~ 1 — 2, 
a choice motivated by the optimal detection window of 
planned laser interferometers such as eLISA. The choice 
of galaxy types and sizes, as well as of the masses of the 
SMBHs just reported, is thus made based on the char- 
acteristic densities expected in the concordance ACDM 
model for galaxies having the same circular velocity as 
the Milky Way V c ~ 200 km s _1 . 

In this study we ch oose a snapshot of th e system at 
the time 2.56 Gyr of ICallegari et al.l (|2011l ). when the 
separation between the two black holes is 700 pc and 
the black holes are still embedded in the individual cusps 
(Figured]). The masses of SMBHs are M.i = 1.5xl0 6 M© 
and M,2 = 4.9 x 10 5 M© corresponding to a mass ratio 
q = M.2/M.1 = 0.3. We select all particles in the central 
5 kpc region of the merger remnant and added all those 
particles, which have pericenter passage smaller than 3 
kpc. 

Most of the gas in the central region is already con- 
verted into stars. The remaining gas particles, which con- 
tribute only a few percent in mass in the central region 
(Figure H]) are also treated as stars in our simulations. 
The total mass of the selected sample is 3.3 x 10 10 M Q . 

Although state-of-the-art hydrodynamical simulations, 
the Callegari et al. (2011) mergers use a few million par- 
ticles to represent the stellar component of a galaxy, it 
is still several orders of magnitude smaller than the ac- 
tual number of stars in a real galaxy. Because of the 
small number of dark matter particles the mass of one 
particle in the primary galaxy is compar able to the mass 
of the SMBH in the satellite galaxy in Callegar i et al.l 
(|2011|) . This is a potential problem for continuing the 
calculation to higher resolution in a regime in which the 
interaction with the background of stars and dark mat- 
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Fig. 2. — Cumulative mass profiles of the different components 
centered at the S MBH M.i of the more massive galaxy at time 
2.56 Gyr (Callcsari et al. 2011). This corresponds to the same time 
when we select our particle sample for direct TV-body simulations 
(T = 0). 

ter might dominate, since it can lead to unrealistic mass 
segregation of the dark matter and non-physically large 
kicks of the two SMBHs. To avoid such high mass parti- 
cle encounters with the SMBH binary, we split each dark 
matter particle in the primary galaxy into ten particles. 
Each child particle has mass 1/10 of the parent particle. 
The child particles are randomly distributed over a 10 
pc sphere corresponding to the softening of the parent 
dark matter particle used in our simulations (see Sec- 
tion 3.2.). The child particles have the same velocities 
as their parent particles. A similar technique of particle 
sp litting, which conser ves linear momentum, was applied 
bv lMaver et~aTl (j2010f ) for SPH particles. We tested that 
split particles do not lead to artificial clumping. 

We change the units of the original model from physical 
units to some model units. We choose our model units, so 
that the total mass of the galaxy (3.3 x 10 10 M Q ) M ga i = 
1, the length unit R is 1 kpc and also G = 1. This results 
in a time unit of 2.6 Myr, velocity unit of 376 km s _1 
corresponding to a speed of light c = 795 in model units. 

For Run-1, which starts at T = (c orresponding to 
the m erger snapshot at time 2.56 Gyr of ICallegar i et al.l 
(|2011h ). we use N = 3071296 particles and evolve the 
system to the point, where the two cusps are merged 
and a "hard binary" is formed. In order to increase the 
computational speed, we reduce the number of particle 
by selecting those, which have a pericenter passage in- 
side 1 kpc N = 1153984 for Run-2. To resolve the stellar 
encounters with the SMBH binary at pericenter passage 
during the late phase of evolution, we reduce the soft- 
ening of star-black hole interactions (see Section [3]) in 
Run-3. Also for Run-3, we again split all dark matter 
particles having pericenter at less than 50 pc leading to 
N = 1327488. The orbit of the SMBH binary is inte- 
grated with higher accuracy during the late phase of the 
binary evolution in Run-4. 

3. NUMERICAL METHODS 
3.1. Simulation Software 

We use the direct iV-body code ip-GPU with 4 th order 
Hermite integration scheme and hierarchical block time 
steps for our iV-body simulations. 

The code is written from scratch in C++ and is based on 
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an earlie r C version of (p-GR APE[l] designed for GRAPE6a 
clusters (lHarfst et al.l 12007V In the present version of 
the ip-GPV code we use native GPU support and direct 
code access to the CPUs using only the NVIDIA CUDA 
librar^. The multi-GPU support is achieved through 
MPI parallelization. Each MPI process uses only a single 
GPU, but we usually start two or more MPI processes 
per node (to effectively use the multi core CPUs and the 
multi GPUs on our clusters). 

The (/3-GPU code is fully parallelized using the MPI 
library. The MPI parallelization was done in the same 
"j" particle parallelization scheme as in the earlier ip- 
GRAPE code. All the particles are divided equally be- 
tween the working nodes (using the MPI_Bcast() com- 
mand) and in each node we calculate only the fractional 
forces for the particles in the current time step, i.e. the 
so called "active" or "i" particles. We get the full forces 
from all the particles acting on the active particles af- 
ter the global MPI_Allreduce() communication routine 
is applied. 

Besides the used 4 th order Hermite integration scheme, 
(p-GP\J additionally supports a 6 th and even 8 th order 
Hermite integration. The numerical integration of the 
particle orbits as well as the time step criterion (see Sec- 
tion 3.3 for det ails) is based on the (seria l CPU) A^-body 
code YEBISU (INitadori fc Makinol (20081) . 

More detail s abou t the ip-GP\J code can be found in 
iBerczik et all (f20lTI ). 

The present version of the code used here is extensively 
modified to handle computational challenges required for 
our current project. 

3.2. Gravitational Softening 

We use a gravitational (Plummer-) softening between 
all particles. ip-GPV supports the use of different soften- 
ing lengths for different components and even individual 
softening for the particles. The softening between the 
SMBH particles is set equal to zero. But the use of zero 
softening for the stars and dark matter leads to the for- 
mation of tight binaries in the system, which causes an 
enormous slow down because of small time steps required 
to resolve the orbits of th e binaries. ip-GPU does not 
include the regularization (Mik kola fc Aarseth 1998) of 
close encounters or binaries, so we use softening to avoid 
the formation of tight binaries. For the interactions be- 
tween two particles, we adopt the following criteria for 
the gravitational softening: 



4 = ( £ ' + e ')A 



(1) 



where e^h = for SMBHs, e s = 0.01 pc for stars and 
£dm = 10 pc for dark matter particles. 

During the late phase of the SMBH binary evolution 
(Run-3), we reduce the star-BH interaction softening ad- 
ditionally by a factor of ten to resolve the stellar encoun- 
ters during the pericenter passage of the SMBH binary. 

3.3. Time Step Criterion 

Th e applied time step criterion ([Nitadori fc Makinol 
2008) for individual particles is 



1 ftp : //ftp . ari . uni-heidelberg . de/staf f /berczik/phi-GRAPE/ 
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Here p is the order of the integrator, a^ k ' is the k-th 
order derivative of the acceleration and r\ is the time step 
parameter. 

</9-GPU can employ different time step parameters r\ for 
different components (BH, stars and dark matter). We 
adopted r\ = 0.1 for all components. In the late phase 
of the SMBH binary evolution (Run-4), we reduce this 
parameter for SMBHs, r\ = 0.3, to integrate the binary 
orbit with smaller time steps, hence achieving higher ac- 
curacy. 

3.4. GPU Clusters 

The A^-body integrations were carried out on three 
GPU high-performance computing clusters, laohu em- 
ploying 172 GPUs at the Center of Information and Com- 
puting at National Astronomical Observatories, Chinese 
Academy of Sciences, titan having 32 GPUs at the As- 
tronomisches Rechen-Institut in Heidelberg, and accre 
employing 192 GPUs at Vanderbilt University, Nashville, 
TN. 

We used up to 64 GPUs for our runs with a total CPU 
wall-clock time of approximately one year. 

4. SMBH BINARY EVOLUTION 

We start our high resolution run at time (T — 0), when 
the two SMBHs are still embedded in their respective 
galaxy cusps. Figure Q] shows the dark matter (top) and 
stellar volume mass densities (bottom) with view on the 
x — y (left) and y — z (right) plane. For visualization 
we use the open-source glnemo2 software packagfl Two 
high density regions around the two SMBHs are clearly 
visible in the figure which suggest that individual galaxy 
cusps are still in the process of merging at the beginning 
of our run. The evolution of the SMBH binary can be de- 
scribed in the following three distinct phases as discussed 
in the following subsections. 

4.1. Dynamical Friction 

In the first phase, the two black holes centered in their 
respective galaxy cusps are unbound to each other and 
move independently in the potential of the merger rem- 
nant. Dynamical friction against the background dark 
matter and stars is very efficient in bringing the two 
SMBHs closer. At about T = 40 Myr the individual 
cusps are already merged into one and the two SMBHs 
are located in a single cusp (see Figure |3J) . Earlier stud- 
ies show that SMBHs form a binary when their relative 
separation ARbh reaches ~ rh (th is the gravitational 
influence radius defined as the radius of a sphere around 
the two black holes enclosing stellar mass equal to twice 
the SMBHs masses). Figure 0] shows the evolution of 
the binary separation. At the same time when the two 

3 http : //pro jets . oamp . f r/proj ects/glnemo2 
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Fig. 3. — Same as in Figure 1 but at T ; 
are embedded in a single cusp. 
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Fig. 4. — Relative separation between the two SMBHs as a 
function of time. The red arrow shows the estimated value of the 
influence radius r^. 



cusps merge (T = 40 Myr) the separation between the 
two SMBHs is roughly about rh — 15 pc and a SMBH 
binary system is formed. 

4.2. Stellar-Dynamical Hardening 

The subsequent evolution of the binary is governed 
by stellar encounters, i.e., predominantly three body en- 
counters. For spherical galaxy models the subsequent 
hard ening is reported to depend on the particle num- 
ber |MiGnoXSnat3 [HI EerczIEeEal] E005). For 
a realistic particle number N, which is several orders 
of magnitude larger than current state-of-the-art simula- 
tions can accommodate, the binary should stall when its 
semi-major axis a ~ a^ for these models. Here ah is the 
semi-major axis of a "hard binary" , as defined by 



a h 



'h 



(1 + a) 2 4 



(4) 



(e. g. lMerritt et~aT1 12007ft . 

In our model r^ is about 15 pc, ah ~ 0.66 pc with 
a^ 1 ~ 1.5 pc -1 . 

Figure [5] shows the evolution of the SMBH binary's 
inverse semi-major axis and eccentricity. The eccentric- 
ity grows to a very high value of e « 0.9 soon after 
the formation of the SMBH binary. The inverse semi- 
major axis of the binary evolution is roughly constant in 
time and the phase of "hard binary" evolution is reached 
quickly. Th i s beha vior is consistent with the findings of 
iKhan et al.1 (J2012I ). who noticed that for shallow cusps 
(with 7 = 0.5, 1.0) l/a of the binary evolves at a con- 
st ajrt_^te_2n™ediately after its formation (see Figure 2 
of lKhan et al.l (|2012f) ). In fact, we find for stellar density 
profile, 7 = 1 = const in our simulation. For steep cusps 
(7 = 1.5, 1.75) the SMBH binary undergoes a rapid phase 
of evolution before entering the hard binary regime, pre- 
sumably corresponding to the clearing of the loss cone 
(|Yull2002t ). The long term evolution of the SMBH binary 
is discussed later. Here we describe in detail the different 
runs that we carried out in order to reach a SMBH sepa- 
ration where gravitational wave emission starts becoming 
important. 

At T = 200 Myr, we reduce the particle number from 
~ 3 million to ~ 1.15 million by selecting the particles 
which have their estimated pericenter in the inner 1 kpc 
to increase the computational speed (Run-2). In order 
to see whether or not our selection has introduced some 
changes in the mass distribution, we plot the cumulative 
mass distribution at various time steps after the new se- 
lection of particles (Figure [B]). The cumulative mass pro- 
file looks very stable in the inner parts. Only in the outer 
parts there is a small expansion of the profile as expected 
due to the new cutoff. Also we start the new run (Run- 
2) 20 Myr earlier to see, if the evolution of the binary 
as it happened in the earlier run (Run-1) can be recov- 
ered. Figure [5] shows that both the inverse semi-major 
axis and the eccentricity evolution are well reproduced 
for the period where the two runs overlap. 

We stop Run-2 when the value of the inverse semi- 
major axis is roughly 40 pc -1 or a — 25 mpc. The value 
of eccentricity at the end of Run-2 is e ~ 0.96. For these 
parameters the value of the pericenter r p for the massive 
binary is r p = (1 — e)a = 1 mpc, which is smaller than the 
softening for the star-BH interaction (7 mpc). In order 
to resolve the star-BH encounters accurately at the peri- 
center passage of the binary, we reduced the softening for 
star-BH encounters by an additional factor 10 and start 
a new run (Run-3) at an earlier time of 520 Myr. For 
Run-3, we again split the dark matter particles, this time 
only those having pericenter smaller than 50 pc from the 
center of the massive binary. As for the earlier splitting, 
each dark matter particle is split into ten particles and 
spread over a 10 pc sphere retaining the velocities of the 
parent particles. We employ the new particle splitting to 
avoid the unphysical jumps in the binary semi-major axis 
caused by massive dark matter particles that occur from 
time to time. Again, our new particle splitting does not 
introduce noticeable changes in the central mass profile. 
For Run-4, we reduce the 77 parameter for the SMBHs 
from 7j = 0.1to?7 = 0.1to achieve higher accuracy in the 
integration of the SMBH binary orbit (there are roughly 
10 orbits in one model time unit). We evolve the SMBH 
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binary till about 1.2 Gyrs with Run-4. The inverse semi- 



major axis value is 50 pc 



20 mpc. The ec- 
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Fig. 5. — The evolution of the binary's inverse semi-major axis 
(top) and eccentricity (bottom). The red arrow in the upper panel 
of the figure points to value of 1/a, which corresponds to the esti- 
mated semi-major axis a^ of the hard binary. The thin blue line 
shows the linear fit to estimate the constant SMBH binary hard- 
ening rate. 
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Fig. 6. — Mass profile after the selection of the new particle sam- 
ple at T = 182 Myrs. The black dotted line represents theoretical 
Hernquist model 7 = 1.0. Clearly the cumulative mass profile is 
very stable in the inner kpc. 



ccntricity value is 0.955, which leads to the pericenter 
distance of 0.9 mpc. The binary evolution for Run-2 is 
very similar to Run-3 and Run-4. The binary's inverse 
semi-major axis evolves at a constant rate (top panel 
of Figure O , which is consistent with our earlier stud- 
ies where we followed the evolu tion of SMBH binary b y 
merging two spherical galaxies (Khan ct al. 2011, 2012). 
We fit a straight line to calculate the binary's hardening 
rate s = gj(l/a) in the late phase of Run-3 and Run-4 
(Figure[5]- top panel). The value of the hardening rate is 
115.7 in model units and 44.5 kpc -1 Myr -1 in physical 
units. This value of the hardening rate is similar to those 
obtained by merging two spheric al galaxies (see top panel 
of Figure 8 of lKhan et al.l (|201lD ) having a similar profile 
(7 = 1) as adopted for the gal axy b ulges in the merge r 
study of ICallegari et al.l (l20Tl . In lKhanet all (l20ll . 
the high value of the hardening rate was attributed to 
the non spherical shape of the merger remnant support- 
ing a large fra c tion o f stars on centrophilic orbits (see also 
iBerczik et al.l (|2006f) V The value of the hardening rate is 
approximately 6 times higher for the same N when com- 
pared to the value for similar mass of SMBH s in spher- 
ical ga laxy models (top panel of Figure 3 of iKhan et al.1 
(2011)). For the merger remnant of two late type galax- 
ies under consideration in our current study, we analyze 
the shape by calculating axes ratios defined for a ho- 
mogeneous ellipsoid with same tensor of inertia. Figure 
UJ shows the intermediate to major and minor to ma- 
jor axes ratios at various distances from the center (top 
panel) and also for different times (bottom panel). We 
can see from the top panel of the Figure [JJ that devi- 
ations from spherical symmetry extend all the way to 
the center (few tens of parsec). The merger remnant is 
considerably flattened, which can also be seen from Fig- 
ure [3j when compared to those which result after the 
merger of two spherical galaxy progenitors. The flatten- 
ing increases as we move to larger distances and becomes 
more or less constant about a distance of 1 kpc. We also 
calculate the axes ratio at different times of evolution 
for the merger remnant at a distance of 200 pc from the 
center, as most of the centrophilic orbits are expected to 
come at about this distance. As is seen from the bottom 
panel of Figure [71 the axes ratios remain constant dur- 
ing the whole time of the evolution of the SMBH binary. 
Due to the non-spherical shape of the merger remnant, 
we expect that the SMBH binary should evolve at a con- 
stant rate supported by the centrophilic orbit family of 
stars rather than the relaxation effects alone. Therefore 
it is reasonable to extrapolate our results for the merger 
of late-type galaxies to a realistic number of star parti- 
cles. Hence we can predict the coalescence time for the 
SMBHs using the estimated hardening rates in the stel- 
lar dynamical phase plus those in the GWs dominated 
regime. 

4.3. Relativistic Regime 

At small enough separation whose value depends on 
the mass and eccentricity of the SMBH binary, gravi- 
tational waves extract energy and angular momentum 
efficiently from the binary, thus making its coalescence 
inevitable. 
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Fig. 7. — Top: Intermediate to major and minor to major axes 
ratios as a function of distance from the center of the SMBH bi- 
nary at T = 200 Myr. Bottom: The time evolution of axes ratios 
calculated at a distance of 0.2 kpc from the center of the SMBH 
binary. 



As is shown in earlier studies (jBerentzen et alJ l2009t 
iKhan et"aT1 120121 ) , the estimated coalescence time ob- 
tained using constant hardening rate s in the s tellar dy- 
namical regime and the formula of Peters ( 1964) for hard- 
ening in the gravitational wave dominated regime agree 
remarkably well with T coa i obtained from simulations 
that follow the binary evolution till coalescence using 
post-Newtonian (VAT) terms in the equation of motion 
of the SMBH binary. But these simulations are compu- 
tationally very expensive due to additional VM terms in 
the equation of motion of the binary and small softening 
needed to resolve the star-BH interactions at pericenter 
till the SMBH binary enters in the gravitational wave 
dominated regime. In order to further evolve the binary 
to the full coalescence of the SMBHs, we need to again 
reduce the softening and also add VM terms in the equa- 
tion of motion of the binary, which would increases the 
computational time drastically (by several months) . 

Further evolution of the SMBH binary can be esti- 
mated by 



da 
~dl 



NB 



I da\ 
\ dt /gw 



^ w+ <!> 



GW 



(5) 



The orbit-averaged expressions — including the lowest 
order 2.5 VM dissipative terms — for the rates of change 
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Fig. 8. — The evolution of the binary's inverse semi-major axis 
(top) and eccentricity (bottom). Runl to Run4 describe the evo- 
lution of the SMBH binary during direct TV-body simulations per- 
formed using (/j-GPU. Thin blue line is the numerical solution of 
the coupled equations JSJl to J7]l for the estimate of the SMBH 
binary evolution. 

of a binary's semi-major axis, and ecc entricity due to 
GW emission are given by iPetersI (|1964l ): 
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Our estimates show that already for the parameters of 
the SMBH binary at T ~ 1.2 Gyr, the contribution to 
the hardening of the SMBH binary can be as large as 
10%. This is the reason that we stop our Run-4 at this 
point. 

We now solve the coupled equations §5$ to (|7|) nu- 
merically to follow the SMBH binary evolution. For a 
numerical solution of the coupled equations, the semi- 
major axis of the binary was chosen at a time T ~ 500 
Myr to have a significant overlap with the iV-body evo- 
lution of the massive binary. The eccentricity value was 
chosen to be 0.95 and we assume that the eccentricity 
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remains constant during the stellar dynamical hardening 
phase. This assumption is supported by the eccentricity 
evolution shown in Figure [5] (bottom) which shows that 
the value of e remains more or less constant from time 
T ~ 600 Myr onwards. 

The estimated evolution is shown in Figure [8] We 
can sec that the coalescence time of the SMBH bi- 
nary is T coa i ~ 2.9 Gyrs. The coalesce nce time of 
2.9 G yr , though longer w h en compared to iKhan et al.l 
IJ20H : IPretoet al.l pOll : IKhan et al.l ([20L2J), is still 
short enough to have a handful 1:10 merge r case s for 
the detection with eLISA. From IKhan et al.l (j2012t ). we 
know that binary hardening rates depend strongly on 
the adopted density profile. For steep density cusps 
with an inner power law density index 7 = 1.75, their 
study shows 4 times higher values of s when compared 
to 7 = 1.0. In the current study the adopted density 
profile at the start of the merger simulation was a Hern- 
quist profile, which has 7 = 1. This slope is observed 
in bright elliptical galaxies which host SMBHs having 
masses ~ 10 8 — 10 9 M©. The faint bulges/ellipticals 
which host smaller SMBHs with masses ~ 10 6 - 10 7 M Q 
have typically steep cusps (7 ~ 1.5 — 1.75). It is con- 
ceivable that the prolonged effect of gas dissipation at 
higher mass and spatial resolution in the last stages of 
the merger, beyond the starting point of the direct TV- 
body calculation, would have led to a steeper baryonic 
cusp at small scales (see Mayer et al. 2007 on the depen- 
dence of the inner density profile of merger remnants on 
the numerical resolution of the gas component). In ad- 
dition, the slope of the initial bulge profile in the galaxy 
models could be steeper and yet still consistent with the 
observed distribution of bulge slopes. Hence, we can eas- 
ily expect that typical coalescence times for SMBH bina- 
ri es in gas-rich merg ers can be shorter and comparable 
to lKhanetaII(|2012f ). 

5. SUMMARY & CONCLUSIONS 



Starting from the results of Callcga ri et al.l (|2011[ ). we 
studied the orbital evolution of a pair of SMBHs in a 
minor merger of disk galaxies (with 30% gas fraction in 
the disk), from an initial separation of 60 kpc to a final 
separation of less than a milli parsec (binary's pericenter 
distance). Initially the mass ratio between the galaxies 
and SMBHs is 0.1. During the merger the two SMBHs 
accrete gas and increase their masses in the process. The 
mass of the SMBH in the satellite galaxy increases almost 
8 fold as the gas in the secondary galaxy is funneled 
towards the center due to the tidal force of the primary 
galaxy at each peri-center passage. The perturbations 
produced by the passages of the secondary galaxy are 
not significant for the primary galaxy, so the SMBH in 
the primary galaxy accretes gas steadily and the mass of 
SMBH here grows by a factor of 2. At the end of SPH 
simulations, the mass ratio between the two SMBHs is 
appro ximately q = 0.3 (see Figure 1 of Callcgar i et al.l 

At the start of our direct iV-body simulations, the sep- 
aration between the two SMBHs is roughly 700 pc, and 
the binary has yet to form. Gas particles, which con- 
tribute only a few percent to the mass of the selected 
central region, are treated as star particles. We use par- 
ticle splitting to reduce the mass of dark matter particles 
to avoid both mass segregation and unphysical encoun- 



ters of high mass dark matter particles with the SMBHs. 
Dynamical friction is very efficient in bringing the two 
SMBHs to a separation where they form a binary at a 
separation of roughly 15 pc. The subsequent harden- 
ing, which happens at a constant rate, is governed by 
individual stars interacting with the massive binary. We 
artificially suppress the contribution of dark matter to 
the hardening of the SMBH binary by introducing a 
large softening (ed m = 10 pc). The shape analysis of 
the merger product shows that the system is predom- 
inantly triaxial from the periphery to the center. The 
SMBH binary evolves at a constant rate and the hard- 
ening rate is high, which suggests that the stalling of the 
SMBH binary should not be an issue in realistic galaxy 
mergers such as those considered here. The eccentricity 
is very high (e ~ 0.95) as was observed for the shallow 
density profile (7 < 1.0) galaxy merger simulations per- 
formed in our earlier studies (jKhan et al.ll20lil [2012). 
The dependence on eccentricity of the coalescence time 
under GW emission is T C oai,GW ~ (1 — e 2 ) 7 / 2 . For very 
eccentric SMBH binaries this could easily account for a 
decrease of an order of magnitude. Accordingly it will be 
very important to further investigate the dependence of 
the eccentricity evolution under different values of 7 and 
q. Currently we are carrying out direct iV-body simula- 
tions of galaxy mergers with galaxies having both steep 
and shallow density profiles at the centers together with 
an initial mass function to address this question. 

With our current study we evolve the SMBH binary 
to a separation (0.9 mpc at pericenter of the SMBH bi- 
nary), where the contribution to the hardening rate of 
the SMBH binary due to the emission of GWs becomes 
important (roughly 10%). Using the constant value of 
the har dening rate in the Newtonian regime and the for- 
mula of lPetersl (|1964f ) for GWs emission from two point 
masses orbiting each other, in the relativistic regime, we 
estimate the coalescence time of two SMBHs to be 2.9 
Gyr after the merger of the two galaxies. The coalescence 
time of 2.9 Gyr, although longer when com pared to the 
times obtained for similar mass binaries in IKhan et al.l 
(2012), is still short enough to have a few 1:10 mergers 
of SMBHs in late type galaxy merg ers in the r a nge a t 
which eLISA is most sensitive. From [Khan et al] (|2012f> . 
we know that binary hardening rates depend strongly 
on the adopted density profile. For steep density cusps 
observed at the center of faint bulges/ellipticals, we can 
expect the coalescence times to be much s horter, com- 
parab le to the ones that were obtained in IKhan et al.1 
(2012) for the merger of steep power law density profile 
galaxies. 

The current work should be regarded mainly as proof- 
of-concept, since we have considered only one particu- 
lar initial condition and have not computed directly the 
binary SMBH shrinking in the post-Newtonian phase. 
Nevertheless, it shows for the first time that the coales- 
cence of the two SMBHs on a timescale sufficiently short 
to be astrophysically relevant does indeed take place as 
a result of a quite realistic galaxy merger with previous 
effects of dissipation taken into account. In the future 
we will explore a wider range of initial conditions mo- 
tivated by cosmological simulations of galaxy formation 
and we will carry out the direct computation of the bi- 
nary shrinking to much smaller separations. 
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